This RMarkdown conducts initial time series analyses on two SNOTEL sites in the Chuska Mountains located in the Navajo Nation. It then compares SNOTEL data to remote sensing data such as SNODAS and PRISM to understand which products best account for snow.

Read in packages

First explore NAs

## # A tibble: 7 x 3
##   variable       n_miss pct_miss
##   <chr>           <int>    <dbl>
## 1 snow_depth_mm      44   1.09  
## 2 air_temp_max_C     29   0.718 
## 3 air_temp_av_C      25   0.619 
## 4 air_temp_min_C     25   0.619 
## 5 swe_mm             15   0.371 
## 6 precip_mm           1   0.0248
## 7 Date                0   0

ENSO

## Warning: Removed 1 rows containing missing values (geom_path).

blue = la nina years

red = el nino years

## Picking joint bandwidth of 25.1

Compare to SNODAS data (extracted for the same lat/long as snotel site)

Daily data

## Joining, by = "Date"
## Picking joint bandwidth of 19.6

## Picking joint bandwidth of 69.3
## Warning: Removed 5 rows containing non-finite values (stat_density_ridges).

## Warning: Removed 2049 rows containing missing values (geom_path).

## Adding missing grouping variables: `month`

## Warning: Removed 2055 rows containing missing values (geom_path).

## Warning: Removed 5 rows containing non-finite values (stat_smooth).
## Warning: Removed 5 rows containing missing values (geom_point).

## Warning: Removed 5 rows containing non-finite values (stat_smooth).

## Warning: Removed 5 rows containing missing values (geom_point).

Timing of max swe

Monthly max from Nov-Apr

## Warning: Removed 3 rows containing non-finite values (stat_smooth).
## Warning: Removed 3 rows containing missing values (geom_point).

Comparison with PRISM data

## Parsed with column specification:
## cols(
##   Date = col_date(format = ""),
##   ppt = col_double()
## )
## Parsed with column specification:
## cols(
##   Date = col_date(format = ""),
##   tmean = col_double()
## )

freezing days

## Picking joint bandwidth of 20.6

Temperature of hottest average temperature days each water year

Number of freezing days, when average temperature is at or below 0

freezing days correlation to SNOTEL swe max

## 
##  Pearson's product-moment correlation
## 
## data:  snotel_swe_max$swe_mm and snotel_swe_max$ndayfr
## t = 3.1496, df = 9, p-value = 0.01174
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2196123 0.9230347
## sample estimates:
##       cor 
## 0.7240951
## 
## Call:
## lm(formula = snotel_swe_max$swe_mm ~ snotel_swe_max$ndayfr)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -110.40  -75.27  -11.91   37.63  187.91 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)  
## (Intercept)           -231.049    151.128  -1.529   0.1607  
## snotel_swe_max$ndayfr    6.491      2.061   3.150   0.0117 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 107.3 on 9 degrees of freedom
## Multiple R-squared:  0.5243, Adjusted R-squared:  0.4715 
## F-statistic:  9.92 on 1 and 9 DF,  p-value: 0.01174

Freezing days relationship to week of water year of max swe

Data Median Mean Variance SD
SNOTEL SWE (mm) 45.72 90.31 13273.55 115.21
SNODAS SWE (mm) 23.00 74.12 11970.78 109.41
PRISM Accumulated Snowfall (mm) 60.54 82.59 6924.39 83.21

## 
## Call:
## lm(formula = wc_prism_max_lm$snow_mm ~ wc_prism_max_lm$snotel_swe_mm)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -34.480 -12.115  -2.662   8.722  46.558 
## 
## Coefficients:
##                               Estimate Std. Error t value   Pr(>|t|)    
## (Intercept)                   66.91982   15.68070   4.268    0.00209 ** 
## wc_prism_max_lm$snotel_swe_mm  0.52618    0.05744   9.160 0.00000739 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 26.81 on 9 degrees of freedom
## Multiple R-squared:  0.9031, Adjusted R-squared:  0.8924 
## F-statistic:  83.9 on 1 and 9 DF,  p-value: 0.000007391
## 
## Call:
## lm(formula = wc_snodas_max_lm$snow_mm ~ wc_snodas_max_lm$snotel_swe_mm)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -66.651  -5.405  11.249  14.802  25.752 
## 
## Coefficients:
##                                 Estimate Std. Error t value     Pr(>|t|)
## (Intercept)                    -14.35588   16.60040  -0.865         0.41
## wc_snodas_max_lm$snotel_swe_mm   1.00035    0.06081  16.449 0.0000000505
##                                   
## (Intercept)                       
## wc_snodas_max_lm$snotel_swe_mm ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 28.38 on 9 degrees of freedom
## Multiple R-squared:  0.9678, Adjusted R-squared:  0.9642 
## F-statistic: 270.6 on 1 and 9 DF,  p-value: 0.00000005051
## [1] 29.37284
## Warning: Removed 25 rows containing missing values (geom_point).

## 
##  Pearson's product-moment correlation
## 
## data:  wc_prism_snotel$snotel_av_temp and wc_prism_snotel$tmean
## t = 704, df = 4008, p-value < 0.00000000000000022
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9957247 0.9962218
## sample estimates:
##       cor 
## 0.9959809
## Warning: Removed 1 rows containing missing values (geom_path).

## [1] 67.06937
## 
## Call:
## lm(formula = wc_prism_snotel_wint$prism_snow_mm_accum ~ wc_prism_snotel_wint$snotel_swe_mm)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -62.83 -33.97 -10.17  15.80 196.90 
## 
## Coefficients:
##                                     Estimate Std. Error t value
## (Intercept)                        39.838709   1.522472   26.17
## wc_prism_snotel_wint$snotel_swe_mm  0.581905   0.009901   58.77
##                                               Pr(>|t|)    
## (Intercept)                        <0.0000000000000002 ***
## wc_prism_snotel_wint$snotel_swe_mm <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 46.55 on 1659 degrees of freedom
## Multiple R-squared:  0.6756, Adjusted R-squared:  0.6754 
## F-statistic:  3454 on 1 and 1659 DF,  p-value: < 0.00000000000000022
## [1] 67.06937

## [1] "Adjusted R2 value is"
## [1] 0.8923609

precipitation comparison

Beaver Spring

First explore NAs

## # A tibble: 8 x 3
##   variable       n_miss pct_miss
##   <chr>           <int>    <dbl>
## 1 snow_depth_mm      50    1.28 
## 2 air_temp_min_C     30    0.770
## 3 air_temp_av_C      29    0.744
## 4 air_temp_max_C     29    0.744
## 5 Date                0    0    
## 6 precip_mm           0    0    
## 7 swe_mm              0    0    
## 8 station             0    0

Whiskey creek vs Beaver springs

## Warning: Removed 30 rows containing missing values (geom_path).

timing of peak swe

SNOTEL vs SNODAS

Daily data

## Joining, by = "Date"
## Warning: Removed 2049 rows containing missing values (geom_path).

## Adding missing grouping variables: `month`

## Warning: Removed 2055 rows containing missing values (geom_path).

## Warning: Removed 8 rows containing non-finite values (stat_smooth).
## Warning: Removed 8 rows containing missing values (geom_point).

## Warning: Removed 8 rows containing non-finite values (stat_smooth).

## Warning: Removed 8 rows containing missing values (geom_point).

## 
## Call:
## lm(formula = bs_snodas_snotel_wint$snodas_swe_mm ~ bs_snodas_snotel_wint$snotel_swe_mm)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -195.405  -11.368   -6.668    8.008  260.892 
## 
## Coefficients:
##                                     Estimate Std. Error t value
## (Intercept)                         6.668346   1.112497   5.994
## bs_snodas_snotel_wint$snotel_swe_mm 0.960421   0.007601 126.351
##                                                 Pr(>|t|)    
## (Intercept)                                0.00000000242 ***
## bs_snodas_snotel_wint$snotel_swe_mm < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 39.07 on 1991 degrees of freedom
## Multiple R-squared:  0.8891, Adjusted R-squared:  0.8891 
## F-statistic: 1.596e+04 on 1 and 1991 DF,  p-value: < 0.00000000000000022
## [1] 39.43443

Comparison with PRISM data

## Parsed with column specification:
## cols(
##   Date = col_date(format = ""),
##   ppt = col_double()
## )
## Parsed with column specification:
## cols(
##   Date = col_date(format = ""),
##   tmean = col_double()
## )

freezing days

Temperature of hottest average temperature days each water year

Number of freezing days, when average temperature is at or below 0

# raw prism snowfall
bs_prism_snow <- bs_prism_freeze %>% 
  select(-waterYear) %>% 
  addWaterYear() %>% 
  mutate(month = month(Date)) %>% 
  drop_na(ppt) %>% 
  group_by(waterYear) %>% 
  mutate(prism_snow_mm_accum = cumsum(prism_snow_mm)) %>%  # accumulate prism precip
  mutate(prism_snow_mm_accum = ifelse(month(Date) %in% c(4,5,6,7,8,9,10), 0, 
                                      prism_snow_mm_accum))

freezing days correlation to SNOTEL swe max

## 
##  Pearson's product-moment correlation
## 
## data:  snotel_swe_max$swe_mm and snotel_swe_max$ndayfr
## t = 2.8648, df = 9, p-value = 0.01863
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1549281 0.9124708
## sample estimates:
##       cor 
## 0.6906187
## 
## Call:
## lm(formula = snotel_swe_max$swe_mm ~ snotel_swe_max$ndayfr)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -125.31  -51.99  -11.22   22.42  176.69 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)  
## (Intercept)           -152.328    139.245  -1.094   0.3024  
## snotel_swe_max$ndayfr    5.365      1.873   2.865   0.0186 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 96.15 on 9 degrees of freedom
## Multiple R-squared:  0.477,  Adjusted R-squared:  0.4188 
## F-statistic: 8.207 on 1 and 9 DF,  p-value: 0.01863

PRISM Snow accumulation vs SNOTEL SWE

## 
## Call:
## lm(formula = bs_prism_snotel_wint$prism_snow_mm_accum ~ bs_prism_snotel_wint$snotel_swe_mm)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -109.15  -30.38  -17.05   19.34  248.51 
## 
## Coefficients:
##                                    Estimate Std. Error t value
## (Intercept)                        29.76708    1.95403   15.23
## bs_prism_snotel_wint$snotel_swe_mm  0.79082    0.01279   61.83
##                                               Pr(>|t|)    
## (Intercept)                        <0.0000000000000002 ***
## bs_prism_snotel_wint$snotel_swe_mm <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 54.66 on 1659 degrees of freedom
## Multiple R-squared:  0.6974, Adjusted R-squared:  0.6972 
## F-statistic:  3823 on 1 and 1659 DF,  p-value: < 0.00000000000000022
Data Median Mean Variance SD
SNOTEL SWE (mm) 76.20 101.06 11526.41 107.36
SNODAS SWE (mm) 54.00 93.46 13775.25 117.37
PRISM Accumulated Snowfall (mm) 68.85 98.13 10146.48 100.73

## [1] 63.38916
## 
## Call:
## lm(formula = bs_prism_max_lm$snow_mm ~ bs_prism_max_lm$snotel_swe_mm)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -66.321 -37.817  -8.161  33.669  88.401 
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    75.8424    36.1926   2.096 0.065593 .  
## bs_prism_max_lm$snotel_swe_mm   0.6667     0.1358   4.910 0.000837 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 54.16 on 9 degrees of freedom
## Multiple R-squared:  0.7281, Adjusted R-squared:  0.6979 
## F-statistic:  24.1 on 1 and 9 DF,  p-value: 0.0008365
## 
## Call:
## lm(formula = bs_snodas_max_lm$snow_mm ~ bs_snodas_max_lm$snotel_swe_mm)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -38.226 -27.618  -6.703  30.295  54.183 
## 
## Coefficients:
##                                Estimate Std. Error t value    Pr(>|t|)    
## (Intercept)                    -45.6611    22.5995   -2.02      0.0741 .  
## bs_snodas_max_lm$snotel_swe_mm   1.1732     0.0848   13.84 0.000000227 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 33.82 on 9 degrees of freedom
## Multiple R-squared:  0.9551, Adjusted R-squared:  0.9501 
## F-statistic: 191.4 on 1 and 9 DF,  p-value: 0.0000002272

## [1] 59.23187

## [1] "Adjusted R2 value is"
## [1] 0.6979208

PRISM snow accumulation using SNOTEL temperature data

precipitation comparison

# combine with snotel data
bs_prism_snotel <- bs_prism_precip %>% 
  mutate(Date = ymd(Date)) %>% 
  rename(prism_ppt_mm = ppt) %>% 
  left_join(bs_clean_metric, by = "Date") %>% 
  select(Date, prism_ppt_mm, "snotel_ppt_mm" = precip_mm, snow_depth_mm, swe_mm)